Le but de ce notebook est de vous rappeler toutes les possibilités intéressantes que vous offre Python en TP de Physique-Chimie.
Retour à la table des matières
Dans cette partie, nous allons voir comment résoudre simplement une équation-différentielle :
Revenir au menu de la résolution des équations différentielles
La résolution d'équations différentielles peut se faire efficacement à l'aide de la fonction odeint du module scipy.integrate. Pour l'utiliser, il faut écrire le problème à résoudre sous la forme :
$$ \frac{dY}{dt}=F(Y,t)$$
$$Y(0)=Y_0$$
Si cette résolution semble a priori être uniquement celle d'une équation d'ordre 1, elle permet en fait de résoudre une équation d'ordre $N$ : il suffit alors que $Y$ soit un vecteur à $N$ éléments. Nous l'étudierons dans les exemples.
A chaque fois, il faudra donc:
y=odeint(F,y0,t,args=()) où args correspond aux arguments de $F$ qui ne sont ni $Y$ ni $t$ (et correspondent à des paramètres du problème par exemple). t est le tableau numpy correspondant aux instants où l'on cherche $Y$.Revenir au menu de la résolution des équations différentielles
On cherche à déterminer la tension $s$ aux bornes d'un condensateur lors de la charge d'un circuit $RC$, régi par l'équation différentielle : $$s+\tau \frac{ds}{dt}=E$$ $$s(0)=0$$
Le problème se met donc sous la forme : $$\frac{ds}{dt}=\frac{E-s}{\tau}=F(s,t)$$ $$s(0)=0$$
Définissons $F$ :
def F(s,t,E,tau): # on met les paramètres du problème (E,tau) en dernier
return (E-s)/tau
Exécutons la fonction odeint et traçons le résultat :
from scipy.integrate import odeint
from scipy.integrate import odepack
import numpy as np
import matplotlib.pyplot as plt
s0=0
tau=10e-3
E=5
t=np.linspace(0,10*tau,1000)
y=odeint(F,s0,t,args=(E,tau))
plt.figure()
plt.plot(t,y)
plt.xlabel("Temps en s")
plt.ylabel("Tension s en V")
plt.grid()
plt.title("Charge d'un condensateur")
plt.show()
Revenir au menu de la résolution des équations différentielles
On s'intéresse désormais à l'équation non linéaire d'ordre 2 modélisation les oscillations d'un pendule de longueur $l$, de masse $m$ et soumis à une force de frottement fluide de coefficient $h$ :
$$\frac{d^2 \theta}{dt^2}+\frac{h}{m} \frac{d\theta}{dt}+\frac{g}{l}\sin(\theta)=0$$
On poser les conditions initiales suivantes :
$$\theta(t=0)=\theta_0$$
$$\frac{d\theta}{dt}(t=0)=\omega_0$$
Pour se ramener au problème d'ordre 1 résolu par odeint, on va poser le vecteur :
$$Y= \begin{bmatrix} \theta \\ \dot{\theta}\end{bmatrix}$$
On déduit :
$$\dot{Y}= \begin{bmatrix} \dot{\theta} \\ \ddot{\theta} \end{bmatrix}=\begin{bmatrix} \dot{\theta} \\ -\frac{h}{m} \dot{\theta}-\frac{g}{l}\sin(\theta) \end{bmatrix} = F(Y,t)$$
Définissons $F$ sous Python :
def F(Y,t,h,m,g,l):
return Y[1],-h/m*Y[1]-g/l*sin(Y[0])
Résolvons ensuite l'équation :
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
from math import sin
y0=[3.14/2,8]
m=0.1
l=1
g=9.8
h=0.02
t=np.linspace(0,100,10000)
y=odeint(F,y0,t,args=(h,m,g,l))
# Tracés
plt.figure()
plt.subplot(211)
plt.plot(t,y[:,0])
plt.xlabel("t en s")
plt.ylabel("Angle theta en rad")
plt.title("Révolutions/Oscillations d'un pendule")
plt.grid()
plt.subplot(212)
plt.plot(y[:,0],y[:,1])
plt.xlabel("theta en rad")
plt.ylabel("theta point en rad/s")
plt.title("Portrait de phase")
plt.grid()
plt.show()
On peut maintenant s'amuser un peu plus en animant tout cela : rendez-vous ici
Revenir au menu de la résolution des équations différentielles
Dans le TP sur l'oscillateur quasi-sinusoïdal à Pont de Wien, on montre que la tension de sortie de l'amplificateur opérationnel $s$ était régie par le jeu d'équations suivant :
où :
Voyons ci-dessous comment résoudre ce problème, en posant : $$Y= \begin{bmatrix} s \\ \dot{s}\end{bmatrix}$$ et $$\dot{Y}=F(Y,t)$$
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
#définissons F
def F(Y,t,Q,w0,K,Vsat):
if abs(Y[0])<Vsat:
d2y=-w0**2*Y[0]-w0/Q*(1-K/3)*Y[1]
else :
d2y=-w0**2*Y[0]-w0/Q*Y[1]
return Y[1],d2y
#posons les constantes
Q=1/3
R=1e3
C=1e-7
w0=1/(R*C)
Vsat=13
K=float(input("Donner la valeur de K :"))
#résolvons
t=np.linspace(0,30e-3,1000)
Y0=[0.0001,0.0001]
Y=odeint(F,Y0,t,args=(Q,w0,K,Vsat))
#traçons
plt.figure()
plt.subplot(211)
plt.plot(t,Y[:,0])
plt.grid()
plt.xlabel("Temps en s")
plt.ylabel ("s en V")
plt.title("Démarrage d'un oscillateur à Pont de Wien")
plt.subplot(212)
plt.plot(Y[:,0],Y[:,1])
plt.grid()
plt.xlabel("s en V")
plt.ylabel("ds/ds en V/s")
plt.title("Portrait de phase")
plt.show()
Donner la valeur de K :3.3
Avant d'étudier cette partie, il faut s'assurer d'avoir revu et compris comment résoudre une équation différentielle avec Python. On pourra passer un peu de temps sur l'exemple de résolution du pendule.
C'est bon ? Alors allons-y. Le but va être de fabriquer cette animation :

Choisissons d'abord un angle initial et une vitesse initiale pour le pendule :
th0=float(input("Donner la valeur de l'angle initial en rad :"))
w00=float(input("Donner la valeur de la vitesse initiale en m/s :"))
Donner la valeur de l'angle initial en rad :2 Donner la valeur de la vitesse initiale en m/s :9
Résolvons l'équation différentielle pour trouver $\theta$ et $\dot{\theta}$ à tout instant :
import numpy as np
from scipy.integrate import odeint
from math import sin, cos
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['figure.figsize']=[8.0,8.0]
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import HTML
#définissons la fonction F associée à l'équation différentielle
def F(Y,t,h,m,g,l):
return Y[1],-h/m*Y[1]-g/l*sin(Y[0])
# définissons les paramètres du problème et les conditions initiales
y0=[th0,w00]
m=0.1
l=1
g=9.8
h=0.02
t=np.linspace(0,100,10000)
dt=t[1]-t[0]
#résolvons l'équation différentielle
Y=odeint(F,y0,t,args=(h,m,g,l))
A l'aide de la fonction pos, nous déterminons les coordonnées de l'extrémité du pendule. Nous créeons aussi les listes temps, theta et thetapoint qui grandiront au cours de l'animation et qui permettront de voir les graphes se tracer progressivement.
def pos(theta):
return sin(theta[-1]),-cos(theta[-1])
On peut alors commencer à coder l'animation :
plot, et de les remplir avec set_data à chaque itération).ìnitial initialisant le tracéanimate qui s'exécutera en boucleOn peut ensuite démarrer l'animation. Préparons tout d'abord la figure :
import matplotlib.animation as animation
fig=plt.figure()
fig.clf()
plt.subplot(221)
l01, =plt.plot([],[],'-')
plt.xlim(-2,2)
plt.ylim(-2,2)
plt.subplot(222)
l1, =plt.plot([],[],'.-')#(t,y[:,0])
plt.xlim(0,20)
plt.ylim(np.min(Y[:,0]),np.max(Y[:,0]))
plt.xlabel("t en s")
plt.ylabel("Angle theta en rad")
plt.title("Révolutions/Oscillations d'un pendule")
plt.grid()
plt.subplot(223)
l2, =plt.plot([],[],'.-')#(y[:,0],y[:,1])
plt.xlabel("theta en rad")
plt.ylabel("theta point en rad/s")
plt.xlim(np.min(Y[:,0]),np.max(Y[:,0]))
plt.ylim(np.min(Y[:,1]),np.max(Y[:,1]))
plt.title("Portrait de phase")
plt.grid()
Et réalisons l'animation désormais :
theta=[Y[0,0]]
thetapoint=[Y[0,1]]
temps=[0]
#fonction d'initialisation
def initial():
global theta
global thetapoint
global Y
global temps
a,b=pos(theta)
l01.set_data([0,a],[0,b])
l1.set_data(t,theta)
l2.set_data(theta,thetapoint)
#fonction d'animation
def animate(i):
global theta
global thetapoint
global Y
global temps
global dt
global th0
global w00
temps.append(temps[-1]+5*dt)
theta.append(Y[5*i,0])
thetapoint.append(Y[5*i,1])
a,b=pos(theta)
l01.set_data([0,a],[0,b])
l1.set_data(temps,theta)
l2.set_data(theta,thetapoint)
HTML(animation.FuncAnimation(fig,animate,np.arange(1, int(len(t)/30)),repeat=False,init_func=initial).to_jshtml())
# c'est un peu long alors on s'arrête à 1/30 du temps total de calcul... modifiable si on a le temps...
# attention cela peut être long